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The existence of an excitation gap in the bulk spectrum is one of the most prominent hngerprints 
of topological phases of matter. In this paper, we propose a family of two dimensional Hamiltonians 
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phase and suggest a concrete physical realization by analyzing the effect of magnetic impurities on 
the surface of strong topological insulators. 
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Introduction — The pursuit of new topological phases 
of matter has lead to the discovery of novel quantum 
states and exotic excitations m- The topological clas¬ 
sification of phases relies on the existence of an excita¬ 
tion gap in the bulk spectrum. Nonetheless, exceptions 
exist, where a system may exhibit topological behavior 
also in the absence of such a bulk gap laig. Among 
these exceptions are, e.g., Weyl semi-metal and nodal 
superconductors 013 , which possess topologically stable 
Fermi points or nodal lines respectively. Additionally, it 
was shown in Ref. [8] that the simultaneous existence of 
one-dimensional gapless modes on the edge and gapless 
modes in the two-dimensional bulk may arise in a family 
of two-dimensional Hamiltonians, generated by coupling 
a topological phase to a gapless phase. 

In this paper, we present a scheme of generating an in¬ 
trinsic gapless superconducting phase in symmetry class 
D P [To] with well localized edge states. This phase is 
distinct from the previous proposals. First, it emerges 
from intrinsic degrees of freedom as opposed to relying on 
an engineered coupling between topological phases and 
gapless ones. Second, it may spontaneously form on the 
surface of three dimensional topological insulators, which 
provides a new experimental route for realizing and prob¬ 
ing gapless topological phases. Finally, the chirality of its 
edge modes depends on disorder in an unusual manner. 
We suggest a concrete example of this phase: a Rashba 
two-dimensional electron gas (2DEG) in the presence of a 
modulated magnetization that is proximity coupled to an 
s-wave superconductor. We analyze the spectral and the 
transport properties of this model both for clean systems 
and in the presence of disorder. 

In Ref. [8] , we have found that in the clean case Dirac 
excitations in the bulk coexist with two types of edge 
states, depending on the edge orientation. For most ori¬ 
entations, we have identified “strong” edge states that 
do not hybridize with the bulk states due to energy and 
momentum conservation. Their wave functions remain 
exponentially localized near the edge despite the absence 
of a bulk gap. For some orientations, this is not the case, 
and the edge states wave functions leak into the bulk. 




- 0.1 - 0.05 0 0.05 0.1 

B, 

(b) 


FIG. I: (a) Phase diagram as a function of Zeeman field 
and disorder strength. Phase transitions are shown in 
red and M denotes a metallic phase, (b) Disorder 
averaged bulk thermal conductance of the model in 
Eq. See m for simulation parameters. 


making the edge modes “weak”. Unlike in the previous 
scheme, here, the intrinsic origin of the edge states leads 
to an additional and surprising dependence of the edge 
states chirality on the lattice termination. 

In the phase we analyze, the inclusion of a uniform 
Zeeman field perpendicular to the 2DEG, Bz, leads to 
a gap opening in the bulk spectrum, and the gapped 
system possesses a non-trivial Ghern number whose sign 
depends on the sign of Bz. Thus, the gapless phase at 
Bz = 0 is a transition between two topologically dis¬ 
tinct insulating phases. In the presence of disorder, the 
non-trivial Ghern number of the system implies delocal¬ 
ized edge states and localized bulk states. At the transi¬ 
tion between different Ghern numbers, the bulk gap must 
close, such that the phase diagram in the space of disor¬ 
der and gap-opening perturbation contains a critical line, 
where the bulk states remain delocalized. The topologi¬ 
cal phase transition occurs at Bz = 0 also in the presence 
of weak on-site potential disorder, leading to the phase 
diagram sketched in Eig. Ef- Starting from the gapless 
point in the clean limit and increasing disorder strength, 
the system remains critical; the bulk states remain delo¬ 
calized while the edge states disappear. 

In the following we introduce and study a general 
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scheme that leads to gapless topological superconductors, 
and then provide a concrete physical realization of such 
a system. This realization involves magnetically doping 
the surface of a strong three-dimensional topological in¬ 
sulator (3d TI), and fulfills the necessary requirements 
without a need for fine tuning. 

Model — It was originally shown by Fu and Kane [2] 
that when a region of a 3d TI surface is proximity cou¬ 
pled to a ferromagnet and another neighboring region 
is proximity coupled to an s-wave superconductor, then 
one-dimensional gapless states must exist at the interface 
between these two regions. Both the ferromagnet and 
the superconductor induce a gap in the spectrum of the 
TI surface. The region in space where the gap changes 
its nature, from a magnetic induced gap to a supercon¬ 
ducting induced gap, is the region that hosts the gapless 
mode. Here, we consider a momentum-space-analogous 
scenario in which the nature of the gap changes in the 
two-dimensional Brillouin zone. We show that this con¬ 
struction dictates the existence of gapless excitations in 
the two-dimensional bulk, which are localized in momen¬ 
tum space and extended in real space. 

The scheme we consider is based on a family of two- 
dimensional Hamiltonians of the type, 

H = ( 1 ) 

where: 

1. Ho represents a Hamiltonian of spinful electrons, 
where spin-orbit coupling breaks the degeneracy of 
the two spin directions for a given momentum. For 
concreteness, we take a 2DEG with a Rashba spin- 
orbit coupling 

Hq = [t (2 — cos akx — cos aky) — ji] (To 

-h A {(7x sin akx — o'y sin aky ), (2) 

where t is the hopping amplitude, /i is the chemical 
potential, A is the spin-orbit strength, a is the lat¬ 
tice constant, and the cr^’s are the Pauli matrices 
in spin space. This Hamiltonian has two circular 
Fermi surfaces, an inner and an outer one. 

2. is a Zeeman coupling to a spatially periodic 
magnetization, characterized by a wave-vector Q, 
that opens a gap in parts of the Fermi surface. For 
concreteness, we take, = mcFz cos Qx, with Q = 
2k where kp is the Fermi momentum of the outer 
Fermi surface. 

3. i^sc = A'0k,t'0-k,4, + H.c. is a superconducting s- 
wave pairing, with '0k,cr being the annihilation oper¬ 
ator of a quasi-particle with spin a and momentum 
k, and A is the induced pairing potential. 

In the absence of superconductivity, the periodic mag¬ 
netization, i^rn, defines a new Brillouin zone of size |Q|, 
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FIG. 2: (Golor online) Band structure of Hamiltonian 
0 with (a) periodic boundary conditions, (b) An edge 
along the y direction. Gapless edge states (red, green) 
coexist with the bulk nodes (blue), (c) Typical edge 
state wave functions. 

and opens a gap at some of its edges [12]. For |Q| ^ 2/ci?, 
where kp is the Fermi momentum, an open Fermi-surface 
develops. The effect of superconducting pairing, i^sc, de¬ 
pends on its strength. Strong pairing (|A| > \m\) renders 
the entire Fermi surface superconducting, and destroys 
the effect of the magnetization. Weak superconductiv¬ 
ity, on the other hand (|A| < |m|), leads to the afore¬ 
mentioned situation where one part of the Fermi surface 
is gapped by the magnetization and the other part is 
gapped by superconductivity. Therefore, the nature of 
the gap changes along a path in momentum space that 
follows the original Fermi surface, and the gap is closed 
at the point of change. 

Introducing the Nambu basis, (-^j, the 

full BdG Hamiltonian becomes 

H =Horz + H^to + AaoTx, (3) 

where the r^’s are Pauli matrices in particle-hole space. 
The model obeys a particle-hole symmetry V = cr^Ty/C, 
= 1, with K signifying complex conjugation, and 
therefore belongs to symmetry class D [9], which has a 
topological classification according to Z in two spatial di¬ 
mensions. In the following, we discretize the Hamiltonian 
(|^ on a square lattice of Lx x Ly sites, setting t = a = 1, 
/i = 0 and A = \/2 — 1, which gives Q = f • 

The band structure near the Fermi level for |A| < \m\ 
is shown in Fig. 0- There are two distinct Dirac nodes 
in the spectrum at the edges of the first Brillouin zone. 
Adding a small uniform Zeeman field, leads to a mass 
gap at the two Dirac cones. The gapped system is a class 
D superconductor with a Ghern number C = sigr].{Bz). 
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In contrast, for |A| > \m\ the spectrum is fully gapped 
(not shown), with C = 0 for any \Bz\ < |A|. For = 0 
and in the presence of boundaries, we find weak edge 
states when the boundary is along the x-direction and 
strong edge states for all other orientations of the edge. 
The spectrum of the system with boundaries along the 
^-direction and for A < m appears in Fig. Dd. Beside 
the two bulk nodes at finite ky^ there are well localized 
(strong) chiral edge modes. 

We find that the edge mode properties are not deter¬ 
mined uniquely by the edge’s orientation. The modu¬ 
lated nature of the magnetization additionally leads to 
a dependence on the termination of the lattice. For the 
parameters we chose, the magnetization wave vector is 
Q = (7r/2,0), z.e., the magnetization is periodic along 
the X direction with a periodicity of four sites, and hence, 
there are four different ways to terminate the lattice in 
the X direction. We label them according to the magne¬ 
tization of the last two sites: (0,+), (+,0), (0,—) and 
(—,0). For a different choice of magnetic periodicity, the 
size of the unit cell changes but the physics remains simi¬ 
lar. The dependence of the edge states on the termination 
is shown in Fig. where the spectrum of a system with 
an edge along the y direction is plotted for the four dif¬ 
ferent terminations. Blue points denote bulk states and 
other colors denote the right edge state. The left edge 
state is fixed in the (0, +) termination and is not shown. 
Both the dispersion and the chirality of the edge state 
change as the termination changes, but the chirality re¬ 
mains non-zero, at zero energy, as long as particle-hole 
symmetry is preserved. Notice that for two terminations 
the edge states cross the Fermi level three times, once 
ai ky = 0 and twice, with an opposite velocity, close to 
the bulk nodes. Hence, the chiralty is determined by the 
edge states at the vicinity of the bulk nodes. 

We confirm the dependence of the edge modes chiral¬ 
ity and their contribution to the thermal conductance 
on both edge orientation and lattice termination by per¬ 
forming transport simulations. Using a three-terminal 
geometry, we separate bulk and edge contributions to 
transport, as described in Appendix A. This can be done 
both in the clean limit as well as in the presence of disor¬ 
der, which we model as a spatial variation of the chem¬ 
ical potential in Eq. drawn independently for each 
lattice site from the uniform distribution [—(5, J], with 
6 being the disorder strength. In the presence of weak 
disorder, the chiral edge states contribution is slightly re¬ 
duced due to the hybridization of the edge modes with 
the bulk. For moderate disorder, the edge states at fi¬ 
nite momenta hybridize strongly with the bulk nodes, as 
opposed to the edge states at small momenta. There¬ 
fore, beyond a certain disorder strength, the edge states 
around zero momentum dominate the transport, lead¬ 
ing to a unique inversion of the edge state chirality (see 
Fig. [^) induced by disorder. As the disorder strength is 
further increased, all edge state contributions go to zero 
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FIG. 3: (Color online) (a) The spectrum of a system 
with an edge along the y direction for all possible 
terminations. The colors in the legend denote the right 
edge state, (b) Strong edge contribution, 

Gch = Grl — Glr^ to the thermal conductance as a 
function of disorder strength 6. We use m = 0.1 and 
A = 0.05. For (5 = 0, the chirality of the edge states 
changes depending on the termination and on disorder. 
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FIG. 4: Edge contribution, Gch, to the thermal 
conductance in the absence of disorder ((5 = 0) as a 
function of system size, for the different terminations. 
Here m = 0.1 and A = 0.05. (a) Strong edge states. 
While their chirality depends on the termination, their 
contribution to the conductance approaches unity, (b) 
Weak edge states. Independently of the termination, 
their contribution approaches zero. 


as the system enters a thermal metal phase. 

In the clean case, the contribution of the strong edge 
states to the thermal conductance becomes quantized in 
the thermodynamic limit as Eig. shows. In contrast, 
the contribution of the weak edge states vanishes as a 
function of system size, see Eig. [^. In this model, the 
coupling between the edge and bulk states depends on 
the normal and superconducting parameters of the model 
and therefore cannot be independently tuned. The en¬ 
ergy scales that correspond to these parameters are large 
compared to the topological gap scale, |A| — |m|, and 
therefore the weak edge states hybridize effectively. In 
the presence of a uniform Zeeman field B^, the depen¬ 
dence on the termination disappears and the chirality is 
determined solely by the sign of Bz. In Appendix B, 
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we provide an illustrative model aiming to intuitively ex¬ 
plain the dependance on the lattice termination. We plot 
the bulk thermal conductance of the system as a function 
of Bz and disorder strength S in Fig. Et- The calculated 
phase diagram agrees nicely with the theoretical expec¬ 
tation of Fig. 

Realization — Class D gapped topological super¬ 
conductors may be realized in hybrid semiconductor- 
superconductor structures, by fine-tuning the chemical 
potential and the applied magnetic field [13]. We take 
a different approach, and consider generating complex 
magnetic order on the surface of a 3d TI without the 
need for fine tuning. This may be realized by means of 
magnetic doping on the surface of the topological insu¬ 
lator [T4f[22] . which we model as a lattice of magnetic 
impurities. 

At low energies, the surface of the undoped system has 
a Dirac cone, exhibiting spin-momentum locking. As the 
Fermi level is tuned away from the Dirac point, the Fermi 
surface becomes nearly hexagonal for a broad range of 
chemical potentials [23H26] . a phenomenon called warp¬ 
ing. Using the model of Fu m, the effective Ruderman- 
Kittel-Kasuya-Yosida (RKKY) interaction [28H32j be¬ 
tween the magnetic impurities is 

^RKKY = [ d^Q 5'q 5'^q, (4) 

Jbz 

where h = 1 and Sq is the Fourier transform of the mag¬ 
netic impurity spin, which we treat as classical. Here, 


jX,X' _ 

-'q - 


— B 

2TTVnc. 


A,A' 


E A,. 

^-(q+G) 


— -1 A, 

= -^ 2^X_’ 
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(q+G)’ 


( 5 ) 


where Jo models an isotropic exchange coupling between 
the magnetic moments and the surface excitations, G 
runs over the set of reciprocal lattice vectors of the im¬ 
purity lattice, and Wc is the area of a unit cell. The quan¬ 
tity X is the spin susceptibility of the bare surface without 
magnetic impurities EH [33], which exhibits pronounced 
peaks at the six nesting vectors of the hexagonal Fermi 
surface [34|. From Eq. ^ we find that when the lattice 
constant of the magnetic moments exceeds 27r|/ci?|“^, the 
RKKY interaction strongly depends on the impurity lat¬ 
tice structure and on the chemical potential. This allows 
to engineer the momenta at which the RKKY interac¬ 
tion is peaked, as we exemplify for different lattice struc¬ 
tures in Appendix D. For lattice constants smaller than 
27r\kF\~^ however, the RKKY interaction is not signifi¬ 
cantly altered by the choice of the lattice. In this regime, 
we use a Metropolis algorithm [35] neglecting the con¬ 
tributions to J away from its peaks, and find that the 
ground state magnetization is a spiral wave, as depicted 
in Fig. whose direction and period are determined by 
one of the nesting vectors. 

To estimate the stability against temperature of this 
order, we introduce the spiral order parameter p = 



FIG. 5: (Color online) (a) Setup illustration and the 
spiral spin order in the ground state for an aligned 
hexagonal lattice of the magnetic impurities, (b) 
Stability of the spiral order (p) against the normalized 
inverse temperature both defined in the text. 


4:7tM~^ Si=i where is the number of simulated 
spins and Qi are the nesting vectors. In the spiral ground 
state p = 1, while it vanishes in the thermodynamic limit 
for an inverse temperature /3 ^ 0 [36]. We calculate the 
dependence of p on /3 for M = 16 and find a transition to 
the spirally ordered phase at Pc ~ 2Aa(, defined as the 
value at which p = ll2 (Fig. EE)- Here, C is a material 
parameter characterizing the warping effect, defined in 
Appendix D. To our knowledge, there is no available data 
for the magnetic couplings Jo on the surface of topologi¬ 
cal insulators, therefore, we employ a reasonable estimate 
of Jo ranging from one meV to about one eV. Using these 
values, the typical range of material parameters [231126] . 
and a spacing of the magnetic impurities of roughly Inm 
- which has been achieved and even underbid in recent 
experiments EHISHI we find a transition temperature 
between 0.12K and 30K. 

After proximity coupling the system to an s-wave su¬ 
perconductor, the phase fulfills the necessary ingredients 
introduced in the presentation of the general scheme. For 
additional details regarding the realization, see Appen¬ 
dices C,D, and E. 

Summary — We have provided a general scheme for 
realizing unusual topological superconductors, which si¬ 
multaneously host gapless Dirac modes in the bulk, and 
chiral Majorana edge states for almost all edge orienta¬ 
tions. We analyze the spectral properties and the thermal 
conductance of the system in the presence of disorder and 
small Zeeman fields, and find that the structure of the 
edge states crucially depends on the termination of the 
lattice. Specifically, for certain terminations, disorder in¬ 
duces an inversion of the edge state chirality. Regarding 
a potential realization, we predict magnetically doped 3d 
TI in the warping regime to exhibit the proposed phase 
by magnetic self-organization when proximity coupled to 
an s-wave superconductor. 
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APPENDIX 



EIG. 6: (Golor online) Sketch of the three-terminal 
geometry used in transport simulations. The edge states 
(arrows) give the only chiral contribution to transport. 


In this appendix, we elaborate on the transport sim¬ 
ulations done for the model presented in the main text 
(Sec. A) and provide an illustrative model aiming to intu¬ 
itively explain the edge mode dependence on the lattice 
termination (Sec. B). We give a more detailed discussion 
on the spiral surface magnetization of strong topological 
insulators, including the cutoff calculations of the spin 
susceptibility in Sec. G, the possibility to engineer the 
RKKY interaction in Sec. D, and the stability of the spi¬ 
ral ground state against temperature in Sec. E. 

A. Transport Simulations 

While both gapped and gapless superconductors are 
perfect conductors of charge, their thermal conductance 
is markedly different. Transport calculations are per¬ 
formed by connecting the system to disorder free leads at 
temperatures Tq and Tq + (5T, and computing the scat¬ 
tering matrix. 


when changing from hard-wall to periodic boundary con¬ 
ditions, as was done in Ref. [8]. In this model however, 
different lattice terminations independently change the 
chirality of the edge states, and result in the formation 
of a spurious conducting channel when periodic bound¬ 
ary conditions are applied. Hence, for a typical termina¬ 
tion, the difference in conductance between systems with 
hard-wall and periodic boundary conditions contains con¬ 
tribution from both the edge and the bulk. Therefore, 
this setup can not be used for all the possible termina¬ 
tions. The three-terminal setup of Eig. [^overcomes this 
problem, while being less prone to finite-size effects, since 
boundary conditions are kept fixed. Additionally, unlike 
the periodic boundary conditions technique, it has the 
advantage of modeling an experimentally accessible sce¬ 
nario. 


B. Dependence on Lattice Termination 


^ab = 


'^ab tab 

ff 

^ab ' ab 


(6) 


between any two leads, a and b. This enables us to de¬ 
termine the thermal conductance in the low-temperature, 
linear response regime. Gab = GoTr (tg^^t^^) where Go = 
7r^/cgTo/6h is the quantum of thermal conductance. All 
transport simulations are performed using the Kwant 
code [39] . 

Conventional topological superconductors have a 
gapped bulk, such that the thermal conductance is only 
due to edge state transport. In contrast, the model intro¬ 
duced in the main text has both bulk and edge excitations 
at the Eermi level, so both contribute to the conductance. 
In order to separate the bulk and edge contributions, we 
perform transport simulations in a three-terminal setup, 
as shown in Eig. By subtracting the conductance be¬ 
tween any two leads from that in the reverse direction 
{Glr — Grl for instance), we obtain the chiral contribu¬ 
tion to transport, which in this model is only due to the 
edge states. 

Transport through the bulk and edge may also be 
determined by comparing the conductance of a system 


The dependence of the strong edge states on the mag¬ 
netic termination may be understood intuitively from the 
small |Q| limit. It was shown by Sau et al US], that the 
combination of a Rashba spin-orbit coupling, s-wave su¬ 
perconducting pairing. A, and a uniform Zeeman field, 
Bz , leads to a two-dimensional class D topological super¬ 
conductor with a well defined Ghern number. The value 
of the Ghern number is sign(H; 2 ) for |A| < \Bz\ and zero 
for |A| > \Bz\. In the limit of small |Q|, the Hamiltonian 
of Eqs. (1-3) in the main text can be locally thought of 
as a system subjected to a uniform Zeeman field. Thus, 
we can think about the system in real space as stripes of 
superconductors, with alternating Ghern numbers, con¬ 
nected in parallel. Edge states flow where the Ghern 
numbers change. Ignoring the gapped bulk states of 
each stripe, we may view the system as a collection of 
one dimensional chiral Major ana channels (see an illus¬ 
tration in Eig. [^. We assume only a nearest neighbors 
coupling between these chiral channels, such that chan¬ 
nels with the same chirality are coupled by tg (symmetric 
coupling) and that channels with an opposite chirality 
are coupled by G (anti-symmetric coupling), and obtain 
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C = 1 C = -l C^l 





It 






FIG. 7: Illustration of the connection of our findings to 
Sau’s model m in the limit of small Q. Upper part: 
the system, in real space, is described as stripes of 
superconductors, with alternating Chern numbers, 
connected in parallel. The chiral Majorana edges states 
of each stripe are also depicted. Lower part: the gapped 
bulk states of each stripe are ignored , and the system 
may be thought of as a collection of one dimensional 
chiral channels connected in parallel and coupled by 
nearest neighbor coupling. 


the two dimensional band-structure of the system. If 
the symmetric coupling is the dominant coupling, i.e., 
ts > ta^ each two neighboring Majorana channels with 
similar chirality form a single Fermionic chiral channel. 
Coupling Fermionic channels with an alternating chiral¬ 
ity through a coupling ta leads to a gapless phase with 
two Dirac cones in the two dimensional Brillouin zone. 
Then, the gapless bulk states in our original model can be 
understood as the band structure emerging from tg > ta. 
The strong edge states are then a remnant of the chan¬ 
nels near the edge, and their nature depend periodically 
on the termination. Consistently with our findings in 
the main part of the manuscript, the chirality of the edge 
states in this picture also depends periodically on the ter¬ 
mination and changes sign with a period of two stripes. 

Adding a uniform Zeeman field on top of this picture 
leads to an asymmetry between regions with different 
Chern numbers. Depending on the sign of the Zeeman 
field, the size of the regions with a given Chern number 
increases while size of the regions with the opposite Chern 
number decreases. In the picture of the one dimensional 
chiral channels, this can be thought of as dimerization of 
channels. Hence, a single Chern number prevails in the 
system. 

It should be pointed out that we provided this cartoon 


model solely as an intuition. Although the model in the 
main text shows similar features to the cartoon model, 
it can not be extrapolated from the the small Q picture, 
since the real space picture in Fig. [^breaks down as Q 
becomes comparable to k-p. 

C. Calculating the Spin Susceptibility with Proper 
Cutoffs 

Following Refs. [211133], we present a detailed calcula¬ 
tion of the surface spin susceptibility in three-dimensional 
topological insulators. We pay special attention to the 
choice of the energy cutoff. 

Model 

The two dimensional surface of a three-dimensional 
topological insulator is well described by the low energy 
Hamiltonian m 

l-Lo = J (fkcl{vo{kx(Ty - kyax) ^jw{k)az)c\^, (7) 

where (Tx,y,z are the Pauli matrices in spin space and 
u;(k) = {k^ + k^)/2 with k± = kx ^ iky. Here, vq is 
the electron velocity near the Dirac point, originating 
from Rashba spin-orbit coupling, and 7 is the warping 
parameter due the cubic Dresselhaus spin-orbit coupling 
of the bulk. We choose the basis Ck = (ck,t,Ck 7 )^ with 
Ck,cr a fermionic annihilation operator for excitations with 
momentum k and spin a. 

In the following, we set the unit of energy to Eq = 
the unit of momentum to ko = and 

h = 1. Also, for simplicity, we express all energy scales 
relative to the energy of the Dirac point. As a function 
of chemical potential, the circular Fermi surface near the 
Dirac point becomes hexagonal, and then later develops 
a ’snowflake’ shape. This phenomenon is called warping 
m, as shown in Fig.[^ For chemical potentials between 
0.55F^o and O.QF^o, marked by red contours in Fig. the 
Fermi surface is almost hexagonal [33] . 

In the following, we place a lattice of magnetic im¬ 
purities on top of the surface, which may be achieved 
experimentally by means of atomic force microscopy. 
Each lattice position is numbered by two integer indices, 
(ji, ^ 2 ) = j ^ The spin of the impurity at rj is de¬ 
scribed by the operator S^. Magnetic moments couple to 
the spin density of the itinerant excitations, = cjcr^cj, 
by a local exchange interaction at position j, 

'Hint = ^ ^ (^) 

xe{x,y,z} 

where are the exchange coupling constants. For sim¬ 
plicity, we assume that all moments couple equally to 
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FIG. 8 : (Color online) Fermi surfaces of Ho in Eq. Q 
for various chemical potentials. The initially circular 
Fermi surface becomes hexagonal with the nesting 
vectors The hexagonal regime is marked by 

red lines. 


the spin density, and that J is spatially homogeneous: 


RKKY interaction 


If the timescales associated to electron dynamics are 
considerably shorter than those of the impurities, the 
electron spin-density operator can be approximately re¬ 
lated to the instantaneous spin configuration of the im¬ 
purities via linear response theory. Then the system is 
well approximated by an RKKY Hamiltonian of the form 

'HrkKY = Yj Y 

X,X'^{x,y,z} 


Here, the non-interacting susceptibility y may be ex¬ 
pressed in terms of the eigenenergies (± 6 k) and eigen¬ 
states of Eq. 0, 




/ oo 

d\ 
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AY. 

-OO ^ ^_ 
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T,p=:k 


T6k p6k-|-q 




( 10 ) 


( 11 ) 


where Ck = ^Jw{kY + is the Fermi-Dirac distri- 

bution function at an inverse temperature (3 with chem¬ 
ical potential {jl^ and 


YE'ik,k') 

= (d.k'^^V’p,k')(d,k'^^Vr,k')*- 

( 12 ) 

We define 



> 

1 

> 

+ 

II 

1 if Te{k),pe{k + q) G (A_, A_|_) 

0 otherwise 

(13) 


in terms of energy cutoffs A±, ensuring the validity of the 
low energy Hamiltonian 0 . 


Evaluating the spin susceptibility 


The integral in Eq. (11) contains a sum over four terms. 


which we label according to the summation indices r and 

p as ++, ±, =F, and-. The ++ contribution originates 

from properties close to the Eermi energy; the ± and ^ 
contributions describe high energy processes between the 

upper and lower branches of the Dirac cone; and the- 

contribution describes processes of states that he deep 
in the Eermi sea. As such, the latter does not play an 
important role at low temperatures. Additionally, since 
bulk states are not close to the Eermi energy nor are 
they localized at the surface, we do not expect them to 
contribute significantly to the RKKY interactions. 

If the RKKY interaction is dominated by the proper¬ 
ties of the system close to the Eermi energy, the ± and ^ 
terms may be neglected, as done in Refs. [miss]. This 
leads to the spin susceptibility shown in Eig. where 
only the ++ term is considered. The dominant contri¬ 
butions occur close to the six nesting vectors (Eig. [^. 
The data is obtained by numerically integrating Eq. ( [ll] ) 
over a 1600 x 1600 hexagonal lattice with a resolution of 
100 X lOO/Zcg. The integrand of Eq. (10) is made well- 
defined at the boundaries of the integration region by 
regularizing the denominator with 77 = 2~^^Eq. 

The mathematical structure of the ±/t contributions, 
however, does not allow to neglected them a priori. In 
fact, without the introduction of the cutoffs A±, the 
spin susceptibility would be dominated by these contri¬ 
butions. We give numerical estimates for the range of 
cutoffs in which their omission is valid. Eig. shows 
the effect of an increasing cutoff, from 0 to 2 Eo, for 
A = A+ = — A_. The ±/t terms contribute to the sus¬ 
ceptibility at small momenta, such that keeping only the 
contributions close to the nesting vectors is not justified 
for sufficiently large A. 

Recent experimental data [231126] shows that for most 
materials the upper cutoff is larger than the absolute 
value of the lower one and that both lie within the 
range where the ±/=F contributions may be omitted. In 
Bi 2 Te 2 Se for instance, the lower cutoff is close to zero 
[24l [25]. In GeBi 2 Te 4 however, both cutoffs are larger 
than for most materials [24]. 


D. Engineering Peak Positions in the RKKY 
Interaction 

The position of the peaks that dominate the RKKY 
interaction in momentum space can be engineered by 
choosing different lattice structures, lattice constants, 
and chemical potentials. The dependence of the RKKY 



















FIG. 9: (Color online) The largest eigenvalue of the 
spin susceptibility y in momentum space at zero 
temperature and for /i = 0.7. 



0.0 0.5 1.0 1.5 

Qxlko 


(a) (b) (c) 

A+ = -A- =0 A+ = -A- = IFo A+ = -A_ = 2Fo 


FIG. 10: (Color online) The largest eigenvalue of the 
spin susceptibility y for different cutoffs A at zero 
temperature and /i = 0.7. For large cutoffs, the 
concentration of y around the nesting momenta gets 
reduced and the low momenta contributions have to be 
taken into account. 


interaction J on the spin susceptibility (Eq. <§ of the 
main text) reads 


jX,X' _ 

- 


— 

27tVuc 


E a,a' 

^-(q+G)’ 

G 


(14) 


where G is the set of the reciprocal lattice vectors and 
Vuc is the unit cell area of the magnetic impurity lattice. 
Given the spin susceptibility, we consider J for different 
lattice structures and lattice constants. For a < 
where Q is the modulus of the nesting vectors, the choice 
of the lattice structure has no significant effect on the 
RKKY interaction. However, when a > 'kQ~^^ differ¬ 
ent peaks overlap since the nesting vectors exceed the 
first Brillouin zone. This allows the position of the peaks 
to be engineered. Exemplary engineered peak positions 
are shown in Eig. There, the largest eigenvalue of 
the RKKY interaction is shown for a a square lattice 
(top panels), a hexagonal lattice with the lattice vectors 
(l,0)a and (l/2,\/3/2)a (middle panels), and a 7r/3 ro- 



-0.4 -0.2 0.0 0.2 0.4 -0.2 -0.1 0.0 0.1 0.2 “O-IO -005 0.00 0.05 0.10 

Hjka q,lko 


(a) a = ^ (b) a = ttQ ^ (c) a = 27 vQ ^ 



qjka 


(d) a = |Q ^ (e) a = ttQ ^ (f) = ^ttQ ^ 



EIG. 11: (Golor online) The largest eigenvalue of the 
RKKY interaction in momentum space for different 
lattice structures and lattice constants, a. The borders 
of the Brillouin zones are marked by dotted lines. A 
square, hexagonal, and 7r/3 rotated hexagonal lattice, is 
shown in the top, middle, and bottom panels, 
respectively. 


tated hexagonal lattice (bottom panels). Changing the 
lattice constant has a similar effect to tuning the chemical 
potential within the hexagonal range. The latter proce¬ 
dure has the advantage of leaving the sample unaltered 
and is more experimentally accessible. 


E. Spiral ground state and Temperature Stability 

In this section, we give details on the Monte Carlo 
algorithm used to determine the transition temperature 
to a spiral ground state. To facilitate the numerics, we 
approximate J by its contributions close the peaks P, 
i.e.. 


jX,X' 

Jp 


E 


rect 





rect 


f % A 
V VA J 


tX,X' 
^ q 5 


(15) 


where A approximates the area of the peak in momen¬ 
tum space and rect is the rectangular function. Erom 
Eig. 1^ we estimate A O.I/cq. A broad class of systems 
with a J-matrix consisting of a finite number of peaks is 
expected to develop a spin wave in their ground states. 
We numerically determine the energetically minimal spin 
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configurations of a 16 x 16 periodic hexagonal lattice of 
spins with a = . Starting from a random spin dis¬ 

tribution where every configuration is equally likely, we 
thermalize the system using a Monte Carlo algorithm [35]. 
For P ^ oo, we reach local energetic minima. We find 
that the initial configurations relax to one of three spiral 
spin waves with momenta Qi, Q 2 , and Qs, cf. Fig. 
One of them is shown in Fig. 5a in the main text. The 
remaining two ground state configurations are related to 
the one in Fig. 5a by a rotation around the 2 : axis about 
an angle of ±7r/3. Note, that each spiral wave has a 
phase degeneracy, equivalent to shifting the origin of the 
spiral. 


Determination of the transition temperature 

The main text introduces the material parameter ( = 
and the inverse transition temperature is given by 

Kq 

Pc ~ 2.4(aC/5c = 2.4 * We estimate values 

of the exchange coupling Jq in order of magnitude ap¬ 
proximation by assigning an exchange energy of the or¬ 
der of 0.1 — leV per unit cell of the topological insula¬ 
tor. This yields J « 10 — lO^meVnm^. We consider 
a lattice constant of a = Inm. Furthermore, we allow 
for k^/EQ ranging from 1.5eV~^nm“^ for Bi 2 Te 2 Se to 
3.8eV~^nm“^ for Bi 2 Te 3 , which we take from the table 
in Ref. [21]. With these values, the critical temperature 
Tc lies between 0.12i^ and 30i^ corresponding to the crit¬ 
ical temperature presented in the main text. 
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